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Abstract 

We study the tracer subgrid term in isopycnal coordinates, Sj. We employ two 

ingredients: the experimental data on vertical spectra of ocean turbulence measured by 

Gargett et al.(1981) and the stochastic approach recently developed by Dukowicz and 

Smith (1997). Our result confirms that Sj is made of two parts: an advection and a 

diffusion term. However, the tracer bolus velocity u consists of two terms 

** 

u = u +u 
1 2 

while in the GM model there is only a term related to u^ which is shown to be: 

u = /cq _1 V a 
1 p 

where q is the thickness weighted average potential vorticity, a result in agreement with 
the recent suggestions by Treguier et al. (1997), Lee et al. (1997) and Greatbatch (1998). 
The second component u^ is new. We compute it in the geostrophic approximation using 
the Gargett et al. data (1981) on ocean vertical turbulence. We find that u >>u t and that 
u^ is orthogonal to u^. 
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I. Introduction 

Ocean general circulation models 0— GCM used in climate studies cannot afford a 
resolution better than a few degrees and therefore are not eddy resolving. This presents a 
difficulty since one must then model the large number of unresolved scales ranging from 
mesoscales down to dissipation scales. Particular attention has recently being paid to the 
mesoscales. Since they are much smaller than their atmospheric analogs, cyclones and 
anticyclones, the assessment of their dynamical role has been considerably more laborious. 
For many years, their role was assumed to be primarily diffusive, a process accounted for 
by a constant and very large diffusivity. But even if one were able to "tune" a set of 
diffusivities to reproduce a particular set of ocean data, such a model would lack predictive 
power and would thus defeat the purpose of being usable in climate modeling where 
prediction of future scenarios is of paramount importance. In addition, and more 
importantly, Gille and Davis (1999) have recently shown that the 'skill index of such 
model is a modest 10% (Z is defined as the percentage of the mean squared mesoscale flux 
reproduced by a given parameterization vis a' vis the value computed with a mesoscale 
resolving model). The advent of eddy resolving models has revealed that mesoscale eddies 
can hardly be represented by a diffusion-type model alone since mesoscale eddies do more 
than smooth out density contrasts, they provide an additional advective transport which, in 
some sense, is the opposite of diffusion since advection brings together regions far removed 
from one another. The advective role of mesoscale eddies was first recognized in 
atmospheric studies due to the work of Andrews and McIntyre (1976), Matsuda (1980), 
Plumb and Mahlman (1987), as reviewed by Andrews et al. (1987). Clearly, a satisfactory 
solution to this problem is possible only within the framework of a complete theory of 
ocean turbulence which entails many interacting scales in a stably stratified medium. But 
even without such a theory, there are observed facts and measured data concerning ocean 
turbulence that may be used to impose strong restrictions on the structure of the subgrid 
term. The most relevant observation is that mixing in the stably stratified ocean occurs 
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preferentially along isopycnal surfaces. This implies that diapycnal fluxes are small and 
thus the non-linear interactions between isopycnal layers are correspondingly weak. 
Neglecting these interactions allows us to conclude that the tracer subgrid terms 
corresponding to isopycnal layers are independent. In other words, they depend only on the 
mean fields characterizing the corresponding layers. Thus, we may view ocean flow as a set 
of two dimensional compressible flows within the layers where the role of "density" is 
played by the thick: less of the layer z^dijdp. Since this feature considerably facilitates 
the solution of the problem of subgrid modeling in isopycnal coordinates, we shall proceed 
to compute the subgrid terms in the iso-coordinates and then to transform them into level 
coordinates, a non trivial process, as we discuss in the following paper. To this theoretical 
argument, one should add the conclusion of Lozier et al.(1994) who showed that averaging 
variables like temperature and salinity over z can give rise to water masses that do not 
exist in the real ocean. 

Progress in modeling subgrids in iso-coordinates has recently been achieved. 
Dukowicz and Smith (1997), Smith (1999) and Dukowicz and Greatbatch (1999) who 
extended to the case of compressible turbulence, the stochastic approach (SA) to 
turbulence developed by Monin and Yaglom (1971). The method can be considered model 
independent as the only assumption is that the diffusion process is Markovian in nature. As 
discussed by Monin and Yaglom (1971), for a diffusion process to be markovian, the time 
considered must be larger than the correlation time of the Lagrangian velocities and if the 
Markovian nature of the process is invalid, the diffusion equation of the Fickian type is 
itself invalid. While the SA is not a closure model in the ordinary sense, it does provide 
considerable insight into the functional dependence of the quantities of interest. In 
particular, the SA showed that the heuristic model of GM (Gent and McWilliams, 1990) 
must be corrected by the presence of new terms which are zero in the GM model but which 
we show below to be non-zero. 

The aim of this paper is take over from where the SA left. We shall supplement it 
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with a new ingredient, the experimental data on ocean turbulence measured by Gargett et 
al. (1981). This will allow us to compute the so-called gauge terms entering the bolus 
velocity which cannot be computed within the SA alone. In this and the next paper, we 
shall show that the new tracer subgrid terms in level and isopycnal coordinates differ from 
those of previous authors. In particular, we recall that in principle one must deal with four 
bolus velocities: 


thickness: u (level), u (iso) 

it ** 

tracer: u (level), u (iso) 

entering the subgrid terms for the mean tracer field r, are (^=level, I=iso): 

= u ++ *V7 : + w' l ’ + r z + Diffusion (level) 


S i = u 


** 


• + Diffusion (iso) 


(la) 

(lb) 

(lc) 

(l d) 


where V(sV H ) is the 2D gradient in level coordinates at constant z, is the 2D gradient at 
constant density p and a tilde/overbar represent averages in iso and level coordinates. In 
previous models, it was assumed that all four bolus velocities are identical 

_]_ * it ** 

u (level)=u (iso)=u (level)=u (iso) (2a) 

whereas we shall show that: 


* ** 
u (iso)=u (iso) 

u + +( level ) = £u"^~ ( level ) 
* 4. 

u (iso)=u (level) + 

where U 2 is a new term orthogonal to and larger than u + . 


(2b) 

(2c) 

(2d) 


II. Tracer Bolus Velocity. Isopycnal Coordinates 

Since the fields under consideration (we omit the x,y,t dependence) 

p{ z), z (P), v(z), v(p) (3a) 

are in general random in nature, we decompose them into a mean and a fluctuating part: 

p( z) = P( z) +p'( z) (3b) 

z (p) = z (p) + z'(p) (3c) 
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(3d) 


For the velocity field, we must use two decompositions, one in level coordinates 

v(z) = v(z) 4- v'(z) 

and the other in isopycnal coordinates: 

v(p) = v(p) + v"(/>) (3e) 

Fluctuations in level and isopycnal are denoted by a prime and a double prime. For z (p) we 
have used the notation (3c) rather than the one corresponding to (3e) since the only 
meaningful average of this function is at fixed p and there is no need to distinguish between 
an overbar and a tilde. In isopycnal coordinates, the transport equation for the tracer 
t- field is written as 

d t r + u-V^r = 0 (3f) 

where 

(3 «> 

where p =dp/d. z. From (3f), we can derive the equation for the mean tracer field t\ 

Li 

^ +s, V +s i =0 (3h) 

where the isopycnal sub-grid term Sj is defined as 

S, s 'u"-y' (3i) 

Since the problem under consideration corresponds to a compressible 2D flow, one should 
consider a "thickness weighted average" defined by 


T = Z ' 1 Z T 

P P 

(4a) 

where z is the thickness defined as: 

p 


z p = 

(4b) 

which, in the case of an adiabatic, Boussinesq fluid, satisfies the equation 


N 

+ 

II 

O 

(4c) 

The average of (4c) gives: 


fa + V^ U) = 0 

(4d) 

tv ^ 

U = u + u 

(4e) 
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u = z'u"/z (4f) 

* 

where u is the ” thickness bolus velocity”. The transport equation for r can be derived 
straightforwardly from Eqs.(3f) and (4c). The result is: 

d^T + u-VpT + Sj = 0 (5a) 

Sj = u* -V^r + R (5b) 


where R is defined by 


R = z'*V • z u"V" (6 

P P P 

where the superscript ("’) indicates fluctuations with respect to the weighted averages, 


T'" = T—T 


Using Eq.(4a), we may write 


t—t = r"— r 1 " = z; 1 z V 1 (7) 

P P 

The SA consists of comparing Eqs.(4d,e), (5) and (6a) with the equations for the time 

evolution of z^ and r which are obtained on the basis of the Fokker-Planck equation and 

which are expressed in terms of the markovian conditional probability density function 

* * 

p(£| £^), where £=(r,t). In Appendix A we derive the following results for R and u : 

R(SA) = -z* l V *(kz V r) - z'*exVf V r (8a) 

v • p p y p p ' p z p p v 

u = Kq-'V^q + 5 p \xV p ip + p X (8b) 


Thus, Sj becomes: 


S = u • V r - z- l V n • (kzVt) 

I p p p v p p > 

* * * 

u = Kq' 1 ? a + f 4 e xV v + z~ l ex\ n tp 
p z p z 


where q is potential vorticity: 


Q = (f+Oz’ 1 


f is the Coriolis parameter, ( is the vertical relative vorticity and e is the unit vector along 

Z 

** 

the z-axis. The velocity u can be interpreted as the " tracer bolus velocity” and in general 

* 

is different from the thickness bolus velocity u because of the presence of 0. The 
diffusivity /t is expressed in terms of the function p(£| £ ) which is considered known within 
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the SA while the functions ip and x arise as arbitrary functions when one equates the 
divergences in Eqs.(4d) and (6a) to the corresponding terms of the SA equations. The 
authors of the SA call ip and \ gauge terms. We find this terminology inadequate since the 
usual meaning from electrodynamcis is that gauge terms are physically arbitrary and no 
physical effect can be attributed to them, which is not the case here. In fact, if we had a 
complete theory of turbulence, the gauge terms could be computed from Eqs.(8a,b) since all 
the other terms would be known. Thus, the functions ip and x ought to be considered 
unknown (since we don't have such a complete theory) rather than arbitrary. Thus, we 
suggest to call \p and \ "axial” terms since they contain vector products of e with a 
gradient which gives rise to an axial vector. Since on the other hand, all terms in (8a-d) 
must be polar vectors, all axial terms must contain the first power of f which is the only 
available pseudoscalar so that the product of an axial vector with a pseudoscalar yields the 
desired polar vector. All other terms in (8a— d) can be called gradient terms. Nevertheless, 
on the basis of the traditional interpretation of gauge terms, Dukowicz and Greatbatch 
(1999) stressed "the rather unsatisfactory situation that our equations depend on gauge" 
and since the tracer subgrid Sj Eq.(8c) depends on x, they suggested the choice: 

X = constant (8f) 

However, although (8f) makes the subgrid Sj independent of the gauge functions, the tracer 
subgrid in level coordinates will still depend on ip. In fact, as shown in paper V, of all the 
functions present in the right hand side of Eqs.(8) only R given by Eq.(8a) is present in 
whereas the other terms of originate from the random nature of the density field. 

HI. Thickness Bolus Velocity. Geostrophic Limit. 

* 

Let us decompose the thickness bolus velocity u (4f) as follows: 

)fC 

u = u + u (9) 

12 v ' 

where 
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(10a) 

(10b) 


u = ( z ' u "y z p = (^z 

U 2 = = “ ^ 

We estimate the two terms as follows: 

|uJ~(z'u")H' 1 (11a) 

|u 2 |~(z T u Tr )(^H)' 1 (11b) 

Here, H is of the order of the ocean depth while <5H is of the order of the correlation length 
scale of vertical turbulence which, following the experimental data of Gargett et al (1981), 
is of the order of few tens of meters (certainly smaller than say 100 meters). Thus, 



£H«H 

(12a) 

and we conclude that 

|u | >> |u I 

(12b) 

* 

1 2 1 


Let us now compute U 2 

in the quasi-geostrophic limit where: 



u" = f^e xV 6" 
z p T 

(13a) 


where <p" is the fluctuating component of the Montgomery potential <f>=gp A Jzdp, and 
satisfies the relation: 


<t>p = gz ' /p ( 13b ) 

Substitute Eqs.(13) into (10b) and take into account that, as we show in paper V, the 
correlation functions of the same fields obtained by averaging in level and isopycnal 
coordinates differ in the second order in the fluctuating fields. This allows us to employ an 
overbar instead of a tilde. The result is as follows: 

u 2 = (2f)- 1 N 2 e z xV p i T2 (14a) 

= (2f)- ! N Vg’VV? (14b) 

= (2f)- 1 NVg- 2 e z x V ? z 2 (14C) 

where N is the Brunt-Vaisala frequency 

N 2 = - gpj p = — glpZp)’ 1 (14d) 

It is worth comparing (14) with the decompositions suggested by Rix and Willebrand 
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(1996, denoted by RW) and by Treguier et al. (1997, denoted by THL). They are: 


u (RW) = (2fN 2 ) _1 g 2 /r 2 e xV/7 2 

2 Z 

(15a) 

u 2 (thl) = f-igV^xV^N- 2 ) 

(15b) 

Using the relation 

P'=~ P/ 

(15c) 

we can further transform (15) to: 

u (RW) = (2fN 2 )‘ 1 e xV (£' 2 N 4 ) 

2 Z 

(16a) 

U 2 (THL) = f- J e z xV (z^N 2 ) 

(16b) 


Both results differ from our expression (14). As for the component of (9), the results of 
THL and RW are identical 


u (THL, RW) = u + = - (?-■ 7TT) Z = (iV) z (16c) 

which are different from our expression (10a). 

IV. Spectra of Ocean Turbulence. The velocity u^ 

Let us rewrite Eq.(13a) in terms of the Fourier components in the isopycnal surface 

which we can approximate by a plane: 

u"(k) = if _1 e xk <^"(k) (17) 

z 

Considering that the gradient of the density is directed almost in the vertical direction and 
that the vertical derivatives of u" and 0" considerably exceed their counterparts in the 
horizontal direction, we obtain: 

oj(k) = if 1 e z xk 0J(k) (18) 

which implies that 

| u "( k )| 2 = rt 2 |^( k )| 2 (w) 

Both horizontal spectra have sharp maxima at 

k =?rL ‘ 1 (20) 

o m v ' 

where L m is the size of the mesoscales. Thus, integrating (19) over k, we obtain the 
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( 21 ) 
(22) 
(23a) 
(23b) 
(24) 

V. The velocity u . Axial terms ip and 

Let us now compare the decomposition (9) with the SA result (8b). We see that u^ has 
the same structure of the last two terms in (8b) and thus we identify them: 

u = z; ] e xV 0 + f _1 e xV y (25) 

It then follows that 

u = /cq 4 V a 
l P H 

u = - + Mf+b-'V^f+b (26) 

This identification has an immediate consequence concerning the choice (8f). In fact, if we 
substitute (25) into (8a) and use (8f), we would obtain a large term in R since the U 2 
obtained from Gargett et al. data (1981) is large because of the small parameter <5H in the 
denominator of Eq.(llb). On the other hand, within an accuracy up to the second order in 
the fluctuating fields, we can rewrite (6a) in the form 

R i/V') (27) 

This contains no derivatives versus z and/or p of the fluctuating fields and thus no 


one-point shear (within a factor of order unity): 


U^ = f 2 k 2 0" 2 
z o z 


The experimental data of Gargett et al. (1981) show that 
As we show in Appendix B, 


N 2 

z 


and thus we adopt the relation: 


u T -u TT | << uT 2 
z z 1 z 


iT 2 = u" 2 
z z 


Substituting Eqs.(21), (22), (23b) and (14d) into (14c), we finally obtain: 

u = ifk* 2 /? e xV i 
2 * o *2 z P P 

The velocity u^ will be discussed in the next section. 


10 



possibility of a <5H which is what makes U 2 large. The choice x=constant would contradict 
(24) and is thus not realizable. 

We can develop the argument further by noticing that when one can neglect the 
derivative of the Coriolis parameter f, both axial functions ip and x can be expressed in 
terms of the only independent variable z . This leads us to conclude that the two terms in 
the right hand side of (25) may differ only by a constant factor and consequently both are 
related to the scale £H. On the other hand, as we have just discussed, ip may be related 
only to the scale H. Thus, the only way to achieve consistency is to take 

ip — constant (28) 

which has the additional implication that the diffusion term R given by (8a) acquires the 
usual form of a diffusion term. Thus, the final expression for Sj is: 

s, = (yy-y - pjp-iwi' y) ( 29 > 

We notice that is orthogonal to and larger than u^. They are given by Eqs.(24) and (26). 
To check that (12b) is satisfied, we take 

MO'V 1 , k* 1 =jr lL m ~ 1 0 4 m, V -L^ 1 (30a) 

where is characteristic horizontal scale. Eq.(24) yields 

u ~10 4 Lj 1 1 (30b) 

A similar estimate of (26) without the /3- term yields (#c~10 3 m 2 s _1 ) 

u^u (GM)~«L^ 1 ~10 3 L^ 1 (30c) 

Thus, the inequality (12b) is fulfilled. 

In summary, since mixing occurs along isopycnals and in view of the fact that the use 
of such coordinates greatly helps the numerical computations (Bleck and Smith, 1990; 
Bleck et al. 1992 and Bleck, private communication), these coordinates are the most 
appropriate. However, since a large majority of ocean models employ level coordinates, it is 
imperative that we study the subgrid in such coordinates, as we do in paper V. 

VI. Comparison with previous work. 
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GM and Gent et al. (1990) have suggested the following heuristic parameterization: 

u*(GM) = -j- , ^(*» / ,i) (31a) 

while Cheniawsky and Holloway (1991, 1993) suggested to consider Eq.(4d) as true 
diffusion of the thickness z p . This implies that rather than (31a) one should use: 

u*(CH) = - B-'y, (31b) 

which has the form of a Fickian diffusion of z p . 

Several authors (McDougall and McIntosh, 1996; Visbeck et al, 1997; Lee et al. 1997; 
THL; Greatbatch 1998; Treguier, 1999; Marshall et al., 1999) have suggested that 
mesoscales exchange potential vorticity rather than thickness and that a model of the form: 

u = «q' ! V q (31c) 

yields better overall results. Eq.(31c) coincides with Eq.(26) for since qsq. When £ is 
unimportant, Eq.(31c) yields 

u = /cf 4 V^f — « z p l V p z p (31d) 

which can be rewritten as 

u* = Kfi A e + u*(CH) (31e) 

where is the unit vector in the y— direction. Eq.(31e) differs from (31b) because of the 
beta— term which becomes important when the gradient of thickness becomes small. For the 

if 

tracer bolus velocity u in Eq.(8c), Greatbatch (1998) suggests the expression: 

u = /cz^V^q (31f) 

which practically coincides with Eqs.(31c-e). On the other hand, if we compare (31d,e) 
with our (26) and (29), we differ by the term u^ Since however the vertical resolution 
required to exhibit such a term is quite high, being contributed by scales of a few tens of 
meters, it is not surprising that its presence may have gone undetected. Much higher 
resolution is required for its detection by eddy resolving models. In this context, it is 
relevant to understand analogous estimates by RW. These authors arrive at a 
decomposition similar to (9) with u coinciding with (16c) while u is expressed by 
Eq.(16a). The discrepancies between our expression (14a) and (16a) is due to the fact that 
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in the latter the random nature of the density function p( z) was not taken into account. 
Even so, a qualitative estimate of (16a) leads once again to (30b). RW used an 
eddy— resolving ocean code, calculated the fluctuating fields in (16a) and (16c) and 
concluded that 

|u| » luj (32) 

while was found to be compatible with GM. An analogous conclusion was reached by 
Treguier (1999). The reason for the discrepancy between (12b) and (32) is that even 
eddy-resolving ocean models are unable to resolve the vertical turbulence that is 
characterized by scales of the order of tens of meters which are the ones that contribute the 
most to (14) and (24). The same reason underlies the discrepancy between our result (24) 
and the result from the linear analysis by Killworth (1997) who considered horizontal scales 
L>>L m and found that the corresponding modes have a vertical length scale ~ H which 
cannot contribute to u^, Eq.(lOb). 

VII. Conclusions 

In this paper we have derived the form of the tracer subgrid in isopycnal coordinates, 

Our main result is Eq.(29): the last term corresponds to the Fickian diffusion while the 

velocities u are given by Eqs.(24) and (26). We note that the u is close to the form (31c) 
1)2 1 

* 

used by Lee et al. (1997) for the whole bolus velocity u . The velocity u 2 is new and it 

corresponds to the axial (gauge) term of the SA results. In addition, u satisfy Eq.(12b). 

1)2 

We have also shown that due to (28), the tracer and thickness bolus velocities coincide. As 
for the subgrid in level coordinates, it must be obtained from the corresponding one in 
isopycnal coordinates. The process is however non trivial since one must take into account 
the random nature of the density field, a problem that will be dealt with in the next paper. 
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Appendix A 


To simplify matters, we adopt a scalar diffusivity 


K if 


(A. 1) 


a simplification that is not essential since at the end we can restore the tensorial nature of 
K. To derive Eqs.(8), we begin with Eq.(15) of Dukowicz and Greatbatch (1999). Using 
(A.l), we have: 


'z = - KZ pV p T + t^xV^r + e z xV p x (A. 2) 

Applying the operator to this equation and taking into account definition (6a), we 
obtain Eq.(8a). Further applying Eq.(A.2) to the case of potential vorticity q and using 
Eqs.(7) and (4f), we transform the left hand side of (A. 2) as follows: 


/ z u'"q" ,V = / z^uq - z^uq = (f+()(u- u ) = _ (f+C) u (A-3) 

Substituting (8e) and (A. 3) into (A. 2) for the case of r=q and using the approximation 
u(=u£ by Dukowicz and Greatbatch (1999), we obtain 

u = «q 4 V p q + z^exS^ + (qz^'^xV^ (A.4) 

where 

X = ~ (x + qtf) ( A - 5 ) 

Changing notation from y-»x in (A.4), we arrive at Eq.(8b). Substituting Eqs.(8a,b) into 
(5b), we arrive at Eqs.(8c,d). 


Appendix B 

To prove (23a), we rewrite (16c) of paper V with an accuracy up to O(p'): 

u'-u" = u ^ - P z \p' (B.l) 

which we further differentiate, square and average. Since £H<<H, it is sufficient to 
differentiate only p'. Thus, we have: 


u'-u" 
z z 


= p~ 2 u 2 p 
y z z y z 


1 2 


(B.2) 


Let us evaluate the derivatives using the length scale H for u z and 6H for p\ We obtains: 
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|u;-u" z |Mh saytpf&p' 2 (b.3) 

Using Eqs. (15c), (13b), (21) and (22), we obtain 

T 2 = ~ (^gA)- 1 ! 2 ^ p~p z (B.4) 

Substituting in (B.3) and taking into account (14d) and (22), we obtain: 

|u;-u" z | 2 / u' 2 = (fLju/jrl^H^H) 2 !? « 10 2 (B.5) 

where we have used Eq.(30a) and the following values: 

N^IO'V 1 , u 2 =10' 3 m 2 s' 2 

<5H 2 =10 3 m 2 , H =10 3 m (B.6) 

Eq.(B.5) proves (23a). 
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